/* 
    extract qcd estimation vs. jet multiplicity in numbers
*/
void qcd_est_by_jets(TH1D *h_qcd_njets, Double_t qcd[5], Double_t qcd_err[5]) {

    int max_njets = sizeof(qcd) / sizeof(Double_t) - 1;

    for ( int i=0; i<5; i++ ) {
        qcd[i]     = 0.;
        qcd_err[i] = 0.;
    }

    TAxis *x = h_qcd_njets->GetXaxis();

    for(int j=1; j <= h_qcd_njets->GetNbinsX(); j++) {

        Int_t   N_jets     = x->GetBinLowEdge(j);
        Double_t N_qcd     = h_qcd_njets->GetBinContent(j);
        Double_t N_qcd_err = h_qcd_njets->GetBinError(j);

        //cout << N_jets << " " << N_qcd << " " << N_qcd_err << endl;

        // 1 jet inclusive
        qcd[0]     += N_qcd;
        qcd_err[0] += TMath::Power(N_qcd_err,2);

        if (N_jets > max_njets) N_jets = max_njets;

        qcd[N_jets]     += N_qcd;
        qcd_err[N_jets] += TMath::Power(N_qcd_err,2); 
    }

    for ( int j=0; j <= max_njets; j++) {
        qcd_err[j] = TMath::Sqrt(qcd_err[j]);
    }

    return;
}

/* 
    recalculate the fake efficiency and its uncertainty
    taking into account the prompt-muon contamination. 
*/
vector<Double_t> eff_fake_contaminate(TString hfile, Int_t njets=0, int mode='+') {

    Double_t ntight = 0.;
    Double_t nloose = 0.;
    vector<Double_t> recal_eff_fake; recal_eff_fake.resize(2);

    TString   tree_eff_fake[4]  = {"data10_eff_fake","mc_ttbar_eff_fake","mc_wjets_eff_fake","mc_zjets_eff_fake"};
    Double_t  scale_xsec_err[4] = {1.0, 0.085, 0.2, 0.2};

    // inverse the scale factor for subtracting mode
    if ( mode == '-' ) {
        for (int i=1; i < sizeof(scale_xsec_err)/sizeof(Double_t); i++) {
            scale_xsec_err[i] = -1.0 * scale_xsec_err[i];
        }
    }

    TFile *f = new TFile(hfile);
    
    for (int i=0; i < sizeof(tree_eff_fake)/sizeof(TString); i++) {
        Double_t num_mu_tight   = 0.;
        Double_t num_mu_loose   = 0.;
        Int_t    num_jets       = 0;
        Double_t d0sig_ctl      = 0.;
        Double_t eff_fake       = 0.;
        Double_t eff_fake_error = 0.;

        TTree *t = (TTree*) f->Get( tree_eff_fake[i] );
 
        t->SetBranchAddress("d0sig_ctl",&d0sig_ctl);
        t->SetBranchAddress("num_jets",&num_jets);
        t->SetBranchAddress("num_mu_tight",&num_mu_tight);
        t->SetBranchAddress("num_mu_loose",&num_mu_loose);
        t->SetBranchAddress("eff_fake",&eff_fake);
        t->SetBranchAddress("eff_fake_error",&eff_fake_error);

        for ( Int_t j=0; j < (Int_t) t->GetEntries(); j++ ) {

            t->GetEntry(j);

            if ( num_jets == njets ) {
                ntight += num_mu_tight * scale_xsec_err[i]; 
                nloose += num_mu_loose * scale_xsec_err[i];
                break;
            }
        }

        cout << tree_eff_fake[i] << " ntight: " << ntight << " nloose: " << nloose << endl;
    }

    // recalculate the fake efficiency
    recal_eff_fake[0] = ( nloose <= 0) ? 0.:( ntight / nloose );

    Double_t dN = nloose - ntight;

    if ( recal_eff_fake[0] != 0. ) {
        recal_eff_fake[1] = TMath::Sqrt(
                                TMath::Power( (nloose + ntight) / TMath::Power(nloose,2) , 2) * ntight +
                                TMath::Power( ntight / TMath::Power(nloose,2), 2) * dN
                            );
    }

    cout << recal_eff_fake[0] << " " << recal_eff_fake[1] << endl;

    f->Close();
    delete f;

    return recal_eff_fake;
}
